Three-dimensional ultrasound imaging but also CT or MR imaging generates three-dimensional data which is usually formatted in a coordinate system which is not a Cartesian coordinate system. Such non-Cartesian grid is e.g. an acoustic grid of an ultrasound probe which obtains images of e.g. the heart of a human being. Such data, which represents a volume, is usually scan-converted along the Cartesian coordinate system before being volume-rendered. When displaying the object scanned, it is necessary to display a surface of the object, wherein such surface shall be displayed “as natural as possible”. One way to do that is to use a virtual external light source and to incorporate the contribution of such external light source into the volume-rendering integral. Such three-dimensional shading is used to render such three-dimensional data set, thereby using application programming interfaces such as OpenGL or DirectX.
For these illumination models such as the Blinn-Phong shading model, it is usually necessary to define a surface normal in order to evaluate the shading model at a specific position. The Blinn-Phong shading model is a modification of the Phong reflection model and is carried out on each point of the surface to be displayed. For this and other shading models it is usually necessary to define the light vector and the surface normal of a specific position of such surface in order to calculate the diffuse reflection and further vectors (such as the eye vector or the camera vector) in order to calculate the specular reflection.
A common way to approximate the surface normal on a specific position is to use a finite differencing scheme such as forward, backward or central differencing. For volumetric data such surface normal is often calculated using the finite difference in all three dimensions (spatial derivatives) also called the spatial gradient vector:
                              g          ⇀                =                  (                                                                                                                f                      ⁡                                              (                                                                              x                            +                            h                                                    ,                          y                          ,                          z                                                )                                                              -                                          f                      ⁡                                              (                                                                              x                            -                            h                                                    ,                          y                          ,                          z                                                )                                                                                                  2                    ⁢                    h                                                                                                                                                                  f                      ⁡                                              (                                                  x                          ,                                                      y                            +                            h                                                    ,                          z                                                )                                                              -                                          f                      ⁡                                              (                                                  x                          ,                                                      y                            -                            h                                                    ,                          z                                                )                                                                                                  2                    ⁢                    h                                                                                                                                                                  f                      ⁡                                              (                                                  x                          ,                          y                          ,                                                      z                            +                            h                                                                          )                                                              -                                          f                      ⁡                                              (                                                  x                          ,                          y                          ,                                                      z                            -                            h                                                                          )                                                                                                  2                    ⁢                    h                                                                                )                                    Formula        ⁢                                  ⁢        1            
The gradient of the scalar field is the direction of highest change, so the surface normal of the isosurface at a certain position is equal to the normalized gradient:
                              n          ⇀                =                              g            ⇀                                                          g              ⇀                                                                      Formula        ⁢                                  ⁢        2            
In order to calculate the gradient of a sampled data volume, a cheap and quick method is to calculate the central differences, thereby using the neighbouring grid values. That is shown in FIG. 1 of the present application:
FIG. 1 shows the apex A of the acoustic probe which is sending out beams 1 such as ultrasound beams in order to obtain two-dimensional cross-sectional images of an object to be scanned, such as the heart of a human being. In order to calculate the gradient G at a certain sample point P (here: only a two-dimensional calculation is shown for the better understanding of the present invention), it is necessary to calculate the two gradient components Ga and Gb, wherein a and b are the axes of a local (finite) coordinate system:Ga=(f(P+x)−f(P−x))/2Gb=(f(P−y)−f(P+y))/2  Formula 3
In a three-dimensional volumetric data set it is clear that it is then necessary to evaluate f, which samples the three-dimensional volumetric data, at six positions. Together with the center sample f(P), which shall be used for the shading task, it is necessary to conduct seven evaluations of f for each shaded position, e.g. each voxel in the 3D voxel space.
It is known in prior art to speed up such calculation by pre-computing the gradient vectors at each discrete position. At each discrete voxel position in a three-dimensional Cartesian coordinate system, the differences are computed and stored together with the original scalar value in a 4-tuple vector (Ga, Gb, Gc, f). In terms of performance, that reduces the number of look-ups when shading a surface within such three-dimensional voxel space by just using the gradient vector components and the scalar value so that it is possible to illuminate these positions by using the pre-calculated 4-tuple vector. The main drawback, however, is the increased data size. The volume data is expanded by at least the factor 4 which easily becomes a problem on machines with a low amount of memory.
For non-Cartesian data, the numerical computation of spatial derivatives (the gradient) and the resultant surface normal becomes more complicated: lighting calculations are based on a global Cartesian coordinate system, wherein light source, viewers position (camera position) etc. are given in such Cartesian coordinate system. Several shading models, such as Blinn-Phong or Phong all use Cartesian coordinate systems for shading calculation.
For non-Cartesian data such as ultrasound data, which is e.g. given in so-called acoustic coordinates, it is not possible to directly use models such as Blinn-Phong or the Phong reflection model thereby using the gradient vector as mentioned above at a given sample position, thereby using any of the differencing schemes mentioned above. This is due to the fact that the coordinate system is not a Cartesian coordinate system. That is also shown in FIG. 1. The local coordinate system with axis a and axis b will not give the correct gradient vector component in direction a, because the values f(P+x) and f(P−x) are not lying in such space. Hence, in order to calculate the gradient vector G, it would be necessary to first re-sample the three-dimensional data set into a Cartesian coordinate system in order to then compute the gradient components, which is then given in the global Cartesian coordinate system.
A different approach is given in U.S. Pat. No. 7,037,263 B2, where the gradient is first determined from data in the acoustic domain and then transformed to the Cartesian coordinate or display screen domain. A transformation matrix is specified for each position within the sampled three-dimensional acoustic data. The above-mentioned solution requires an expensive re-sampling step to a Cartesian geometry.